library(Amelia)
library(Zelig)
library(xtable)
library(emdbook)
library(pscl)
library(fields)

## Load Data ##
# ISQBalcells_corrected.dat: New Dataset containging corrected measure of CNT Affiliation received from author#
balcells.isq <- read.dta("ISQBalcells_corrected.dta")  
# balcells.old <- read.dta("ISQBalcells.dta")
balcells.jcr <- read.dta("CataloniaJCRreplication.dta") # a related publication in JCR 2011

## Data Preprocessing ##
balcells.isq <- balcells.isq[-1063,] # Get rid of irrelevant row
balcells.isq$cntaffiliation <- balcells.isq$cntaffiliation2 ## Replace corrected column
balcells.isq <- balcells.isq[,-(14:15)] # Get rid of unneeded columns

## Create consolidated dataset: balcells ##
jcrvariables <- c("support_left", "support_right", "compabs", "popmil", "bombdum", 
                  "bombkilled10", "latitude", "longitude")
balcells <- cbind(balcells.isq, balcells.jcr[jcrvariables]) 
balcells$latitude[balcells$latitude==0] <- NA ## Change incorrect 0s to NAS ##
balcells$longitude[balcells$longitude==0] <- NA ## Change incorrect 0s to NAS ##

# Remove unneeded objects #
rm(balcells.isq, balcells.jcr,jcrvariables)

## Load User Created Functions ##
source("User-Created Functions.R")

## Set aside clean dataset for multiple imputation (Appendix) ##
impute <- balcells

## Transformation for Data Analysis ##
balcells$support_right <- balcells$support_right/100 ## Rescale Support_Right to 0 to 1 
balcells$support_left <- balcells$support_left/100 ## Rescale Support_Right to 0 to 1
balcells$logpopulation <- log(balcells$population) ## Add logpopulation
balcells$union <- NA
balcells$union[(balcells$cntaffiliation>0|balcells$ugtaffiliation>0)] <- 1
balcells$union[(balcells$cntaffiliation==0&balcells$ugtaffiliation==0)] <- 0
balcells$sr_union <- balcells$support_right*balcells$union

### Subsetting for Data Analysis ###

original <- balcells
union <- balcells[(balcells$cntaffiliation>0|balcells$ugtaffiliation>0),]
nonunion <- balcells[(balcells$cntaffiliation==0|balcells$ugtaffiliation==0),]

#####################################
########## Entire Dataset ###########
#####################################

### I. Full ZINB model, Table 2, Balcells(2010) ###
zinb <- zeroinfl(executed_left ~ competition + frontline + population + cntaffiliation 
                 + ugtaffiliation + border + sea + altitude + catholiccenter | competition 
                 + frontline + population + ugtaffiliation + border + sea + altitude, 
                 data=original, dist = "negbin")

## residuals ##
zinb_original_residuals <- zinb$residuals
dev.off()
plot(zinb_original_residuals)
plot(table(zinb_original_residuals))
stats(zinb_original_residuals)

## Set Values for Simulation ##
data <- original[,3:dim(original)[2]] # cut out the id and names of provinces
data <- apply(data, 2, mean,na.rm=T) # mean everything

count.vars <- c("competition", "frontline", "population", "cntaffiliation", 
                "ugtaffiliation", "border", "sea", "altitude", "catholiccenter") 
zero.vars <-c("competition", "frontline", "population", "ugtaffiliation", "border", 
              "sea", "altitude" )
vary <- "competition"

setx <- XZout_Function(data=data, count.vars=count.vars, zero.vars=zero.vars, vary=vary)

## Simulated results ##
sim.result <- ExpPred_Function.ZINB(zinb, x.out = setx$x.out, z.out = setx$z.out, 
                                    par.simN = 10^5, predict.simN = 10^5)

## Replication of Balcells 2010, Figure 1 ##
complevel <- seq(from = 0, to = 1, by = 0.1)
png("Figure1_Left.png", width=550, height=450)
plot(complevel, apply(sim.result$predict, 2, mean), type="b", ylim=c(1,5.3), xlab="Competition", ylab="Number of Executions")
dev.off()

# Predicted and Expected Value Plots
dev.off()
png("Figure1_Right.png", width=550, height=450)
plot(complevel, apply(sim.result$predict, 2, mean), type="b", ylim=c(0,18), xlab="Competition", ylab="Number of Executions")
#title( main="Figure 1. Predicted Number of Executions by Competition", font.main = 1, cex.main=1)
lines(complevel, apply(sim.result$predict, 2, quantile, .025), type="l", lty="dotted")
lines(complevel, apply(sim.result$predict, 2, quantile, .975), type="l", lty="dotted")

### Expected Values are Similar
#dev.off() 
#plot(complevel, apply(sim.result$expect, 2, mean), type="b", xlab="Competition", ylab="Number of Executions",)
#title( main="Simulated Figure 1 (with 95% Confidence Intervals)", font.main = 1, cex.main=1)
lines(complevel, apply(sim.result$expect, 2, quantile, .025), lty="dashed")
lines(complevel, apply(sim.result$expect, 2, quantile, .975), lty="dashed")
legend(x="topleft", c("Predicted Value CI","Expected Value CI"), lty=c("dotted","dashed"))
dev.off()  

# Range and SD of Executed Left #
stats(balcells$executed_left)

#########################################################
############# Use Logpop with Original Dataset ##########
#########################################################

# Deaths total
plot(density(balcells$executed_left))
dev.off()
png("FigureDeaths.png", width=550, height=450)
plot(table(balcells$executed_left),xlab="Executions",
     ylab="Number of observations")
dev.off()

# Outlier
balcells[balcells$executed_left==max(balcells$executed_left),]

# Deaths per Capita
plot(density(balcells$executed_left/balcells$population*1000,na.rm=T))
dev.off()
png("FigureDeathsPerCap.png", width=550, height=450)
plot(table(round(balcells$executed_left/balcells$population*1000,0)),
     xlab="Executions per 1000 people",
     ylab="Number of observations")
dev.off()

# Outlier
deathspercap <- balcells$executed_left/balcells$population
deathspercap
balcells[deathspercap==max(deathspercap,na.rm=T),]

####### ZINB Model with Logpop*1 #########

### * Note that we use population rather than logpopulation in the zero model,   ###
### to hold everything as similar as possible to her model. Using logpopulation  ###
### does not lead to substantively or statistically signficant changes in results###
zinb <- zeroinfl(executed_left ~ competition + frontline + offset(logpopulation) + cntaffiliation 
                 + ugtaffiliation + border + sea + altitude + catholiccenter | competition 
                 + frontline + population + ugtaffiliation + border + sea + altitude, 
                 data=original, dist = "negbin"); summary(zinb)

## Set Covariates Values for Simulation ##
data <- original[,3:dim(original)[2]] 
data <- apply(data, 2, mean,na.rm=T) 
median.values <- c("frontline","border","sea","catholiccenter") 
data[median.values] <- apply(original[,median.values],2, median,na.rm=T)

count.vars <- c("competition", "frontline", "cntaffiliation", "ugtaffiliation", "border", "sea", "altitude", "catholiccenter") 
zero.vars <-c("competition", "frontline", "population", "ugtaffiliation", "border", "sea", "altitude")
vary <- "competition"
logpop <- mean(original$logpopulation, na.rm=T)[1]

setx <- XZout_Function(data=data, count.vars=count.vars, zero.vars=zero.vars, vary=vary)

## Simulate Results ##
sim.result <- ExpPred_Function.ZINB.log(zinb, x.out=setx$x.out, z.out=setx$z.out, logpop=logpop, par.simN=10^5, predict.simN=10^5)
# Predicted Executions, Support Right #
predict <- apply(sim.result$predict, 2, mean); predict
apply(sim.result$expect, 2, quantile, .025)
apply(sim.result$expect, 2, quantile, .975)
# Expected Executions, Support Right #
expect <- apply(sim.result$expect, 2, mean); expect


### Interaction Right Support & Union / Including Union Variables / Logpop*1 ####
zinb <- zeroinfl(executed_left ~ support_right + union + support_right:union + frontline 
                 + offset(logpopulation) + cntaffiliation + ugtaffiliation + border + sea + altitude + catholiccenter | support_right 
                 + union + support_right:union + frontline + population + ugtaffiliation + border + sea + altitude, 
                 data=original, dist = "negbin") ; summary(zinb)

count.vars <- c("support_right", "union", "frontline", "cntaffiliation", "ugtaffiliation", "border", "sea", "altitude", "catholiccenter", "sr_union") 
zero.vars <-c("support_right", "union", "frontline", "population", "ugtaffiliation", "border", "sea", "altitude", "sr_union")
vary <- "support_right"

## residuals ##
zinb_dummy_residuals <- zinb$residuals
plot(table(zinb_dummy_residuals))
hist(zinb_dummy_residuals)
stats(zinb_dummy_residuals)

# simulated effects
data <- original[,3:dim(original)[2]] 
data <- apply(data, 2, mean,na.rm=T)
median.values <- c("frontline","border","sea","catholiccenter") # choose the median columns
data[median.values] <- apply(original[,median.values],2, median,na.rm=T) # make them median

setx <- XZout_Function(data=data, count.vars=count.vars, zero.vars=zero.vars, vary=vary)

# Union=1 
setx$x.out[,3] <- 1
setx$z.out[,3] <- 1

setx$x.out[,11] <- setx$x.out[,2]
setx$z.out[,10] <- setx$z.out[,2]
logpop <- mean(original$logpopulation, na.rm=T)

sim.result.union <- ExpPred_Function.ZINB.log(zinb, x.out=setx$x.out, z.out=setx$z.out, logpop=logpop, par.simN=10^5, predict.simN=10^5)

# Union=0 
setx$x.out[,3] <- 0
setx$z.out[,3] <- 0

setx$x.out[,11] <- 0
setx$z.out[,10] <- 0
logpop <- mean(original$logpopulation, na.rm=T)
sim.result.nonunion <- ExpPred_Function.ZINB.log(zinb, x.out=setx$x.out, z.out=setx$z.out, logpop=logpop, par.simN=10^5, predict.simN=10^5)

### Figure 2. Expected Number of Executions (Union vs Non-Union) ####
dev.off()
png("Figure2_Left.png", width=550, height=450)
plot(complevel, apply(sim.result.union$expect, 2, mean), lty="solid", type="l", ylim=c(0,20), xlab="Right Support", ylab="Number of  Executions", col=4, lwd=2)
lines(complevel, apply(sim.result.nonunion$expect, 2, mean), type="l", col=2, lwd=2)
legend(x="topleft", c("Union=1","Union=0"), lty=c("solid","solid"), col=c(4,2), lwd=c(2,2))
lines(complevel, apply(sim.result.union$expect, 2, quantile, .025), lty="dotted", col=4, lwd=2)
lines(complevel, apply(sim.result.union$expect, 2, quantile, .975), lty="dotted", col=4, lwd=2)
lines(complevel, apply(sim.result.nonunion$expect, 2, quantile, .025), lty="dotted", col=2, lwd=2)
lines(complevel, apply(sim.result.nonunion$expect, 2, quantile, .975), lty="dotted", col=2, lwd=2)
dev.off()

# Not included Figure. Predicted Number of Executions ######
png("Figure2_Right.png", width=550, height=450)
plot(complevel, apply(sim.result.union$predict, 2, mean), lty="solid", type="l", ylim=c(-5,100), xlab="Right Support", ylab="Predicted Executions", col=4, lwd=2)
lines(complevel, apply(sim.result.nonunion$predict, 2, mean), type="l", col=2, lwd=2)
lines(complevel, apply(sim.result.union$predict, 2, quantile, .025), lty="dotted", col=4, lwd=2)
lines(complevel, apply(sim.result.union$predict, 2, quantile, .975), lty="dotted", col=4, lwd=2)
lines(complevel, apply(sim.result.nonunion$predict, 2, quantile, .025), lty="dotted", col=2, lwd=2)
lines(complevel, apply(sim.result.nonunion$predict, 2, quantile, .975), lty="dotted", col=2, lwd=2)
dev.off()


## Different Specifications : Competition & Union Interaction ##
zinb <- zeroinfl(executed_left ~ competition + union + competition:union + frontline 
                 + offset(logpopulation) + cntaffiliation + ugtaffiliation + border + sea + altitude + catholiccenter | competition 
                 + union + competition:union + frontline + population + ugtaffiliation + border + sea + altitude, 
                 data=original, dist = "negbin") ; summary(zinb)

count.vars <- c("competition", "union", "frontline", "cntaffiliation", "ugtaffiliation", "border", "sea", "altitude", "catholiccenter", "sr_union") 
zero.vars <-c("competition", "union", "frontline", "population", "ugtaffiliation", "border", "sea", "altitude", "sr_union")
vary <- "competition"

data <- original[,3:dim(original)[2]] 
data <- apply(data, 2, mean,na.rm=T)
median.values <- c("frontline","border","sea","catholiccenter") # choose the median columns
data[median.values] <- apply(original[,median.values],2, median,na.rm=T) # make them median

setx <- XZout_Function(data=data, count.vars=count.vars, zero.vars=zero.vars, vary=vary)

# Union=1 
setx$x.out[,3] <- 1
setx$z.out[,3] <- 1

setx$x.out[,11] <- setx$x.out[,2]
setx$z.out[,10] <- setx$z.out[,2]
logpop <- mean(original$logpopulation, na.rm=T)

sim.result.union <- ExpPred_Function.ZINB.log(zinb, x.out=setx$x.out, z.out=setx$z.out, logpop=logpop, par.simN=10^5, predict.simN=10^5)

# Union=0 
setx$x.out[,3] <- 0
setx$z.out[,3] <- 0

setx$x.out[,11] <- 0
setx$z.out[,10] <- 0
logpop <- mean(original$logpopulation, na.rm=T)
sim.result.nonunion <- ExpPred_Function.ZINB.log(zinb, x.out=setx$x.out, z.out=setx$z.out, logpop=logpop, par.simN=10^5, predict.simN=10^5)

### Competition Interaction (not included). Expected Number of Executions (Union vs Non-Union) ####

dev.off()
plot(complevel, apply(sim.result.union$expect, 2, mean), lty="solid", type="l", ylim=c(0,20), xlab="Right Support", ylab="Number of  Executions", col=4, lwd=2)
lines(complevel, apply(sim.result.nonunion$expect, 2, mean), type="l", col=2, lwd=2)
legend(x="topleft", c("Union=1","Union=0"), lty=c("solid","solid"), col=c(4,2), lwd=c(2,2))
lines(complevel, apply(sim.result.union$expect, 2, quantile, .025), lty="dotted", col=4, lwd=2)
lines(complevel, apply(sim.result.union$expect, 2, quantile, .975), lty="dotted", col=4, lwd=2)
lines(complevel, apply(sim.result.nonunion$expect, 2, quantile, .025), lty="dotted", col=2, lwd=2)
lines(complevel, apply(sim.result.nonunion$expect, 2, quantile, .975), lty="dotted", col=2, lwd=2)


dev.off()
plot(complevel, apply(sim.result.union$predict, 2, mean), lty="solid", type="l", ylim=c(-5,100), xlab="Right Support", ylab="Predicted Executions", col=4, lwd=2)
lines(complevel, apply(sim.result.nonunion$predict, 2, mean), type="l", col=2, lwd=2)
lines(complevel, apply(sim.result.union$predict, 2, quantile, .025), lty="dotted", col=4, lwd=2)
lines(complevel, apply(sim.result.union$predict, 2, quantile, .975), lty="dotted", col=4, lwd=2)
lines(complevel, apply(sim.result.nonunion$predict, 2, quantile, .025), lty="dotted", col=2, lwd=2)
lines(complevel, apply(sim.result.nonunion$predict, 2, quantile, .975), lty="dotted", col=2, lwd=2)


##########################################################################
################## UNION / NONUNION SUBSET DATA ##########################
##########################################################################

### Support Right w/ LogPopulation ###
nb <- glm.nb(executed_left ~ support_right + frontline + offset(logpopulation) + cntaffiliation 
            + ugtaffiliation + border + sea + altitude + catholiccenter,
            data=union) ; summary(nb)

## residuals ##
nb_union_residuals <- nb$residuals
plot(nb_union_residuals)
hist(nb_union_residuals)
stats(nb_union_residuals)


x.out <- setx(nb,support_right=seq(0,1,0.1))
x.out$border <- 0
x.out$sea <- 0
x.out$catholiccenter <-0
logpop <- mean(union$logpopulation, na.rm=T)

sim.result <- ExpPred_Function.NB.log(nb, x.out=x.out, logpop=logpop, par.simN=10^5, predict.simN=10^5)

png("Figure3.png", width=550, height=450)
plot(complevel, apply(sim.result$predict, 2, mean), type="b", ylim=c(0,100), xlab="Right Support", ylab="Number of Executions")
#title( main="Figure 1. Predicted Number of Executions by Competition", font.main = 1, cex.main=1)
lines(complevel, apply(sim.result$predict, 2, quantile, .025), type="l", lty="dotted")
lines(complevel, apply(sim.result$predict, 2, quantile, .975), type="l", lty="dotted")

lines(complevel, apply(sim.result$expect, 2, quantile, .025), lty="dashed")
lines(complevel, apply(sim.result$expect, 2, quantile, .975), lty="dashed")
legend(x="topleft", c("Predicted Value CI","Expected Value CI"), lty=c("dotted","dashed"))
dev.off()  

### Different Specifications: Competition w/ LogPopulation ###
nb <- glm.nb(executed_left ~ competition + frontline + offset(logpopulation) + cntaffiliation 
             + ugtaffiliation + border + sea + altitude + catholiccenter,
             data=union) ; summary(nb) ### Competition not statistically significant

x.out <- setx(nb,competition=seq(0,1,0.1))
x.out$border <- 0
x.out$sea <- 0
x.out$catholiccenter <-0
logpop <- mean(union$logpopulation, na.rm=T)

sim.result <- ExpPred_Function.NB.log(nb, x.out=x.out, logpop=logpop, par.simN=10^5, predict.simN=10^5)

plot(complevel, apply(sim.result$predict, 2, mean), type="b", ylim=c(0,60), xlab="Competition", ylab="Number of Executions")
#title( main="Figure 1. Predicted Number of Executions by Competition", font.main = 1, cex.main=1)
lines(complevel, apply(sim.result$predict, 2, quantile, .025), type="l", lty="dotted")
lines(complevel, apply(sim.result$predict, 2, quantile, .975), type="l", lty="dotted")

lines(complevel, apply(sim.result$expect, 2, quantile, .025), lty="dashed")
lines(complevel, apply(sim.result$expect, 2, quantile, .975), lty="dashed")
legend(x="topleft", c("Predicted Value CI","Expected Value CI"), lty=c("dotted","dashed"))

#### Nonunion Subset Data #########

### Competition & Logpopulation ####
zinb <- zeroinfl(executed_left ~ competition + frontline + offset(logpopulation) 
                 + cntaffiliation + ugtaffiliation + border + sea + altitude  
                 + catholiccenter | competition + frontline + population + ugtaffiliation + border + sea + altitude, 
                 data=nonunion, dist ="negbin"); summary(zinb)

data <- nonunion[,3:dim(nonunion)[2]] 
data <- apply(data, 2, mean,na.rm=T) 
median.values <- c("frontline","border","sea","catholiccenter") 
data[median.values] <- apply(nonunion[,median.values],2, median,na.rm=T)

count.vars <- c("competition", "frontline", "cntaffiliation", "ugtaffiliation", "border", "sea", "altitude", "catholiccenter") 
zero.vars <-c("competition", "frontline", "population", "ugtaffiliation", "border", "sea", "altitude")
vary <- "competition"
logpop <- mean(nonunion$logpopulation, na.rm=T)[1]

setx <- XZout_Function(data=data, count.vars=count.vars, zero.vars=zero.vars, vary=vary)

sim.result <- ExpPred_Function.ZINB.log(zinb, x.out=setx$x.out, z.out=setx$z.out, logpop=logpop, par.simN=10^5, predict.simN=10^5)

plot(complevel, apply(sim.result$predict, 2, mean), type="b", ylim=c(0,15), xlab="Right Support", ylab="Number of Executions")
#title( main="Figure 1. Predicted Number of Executions by Competition", font.main = 1, cex.main=1)
lines(complevel, apply(sim.result$predict, 2, quantile, .025), type="l", lty="dotted")
lines(complevel, apply(sim.result$predict, 2, quantile, .975), type="l", lty="dotted")

lines(complevel, apply(sim.result$expect, 2, quantile, .025), lty="dashed")
lines(complevel, apply(sim.result$expect, 2, quantile, .975), lty="dashed")
legend(x="topleft", c("Predicted Value CI","Expected Value CI"), lty=c("dotted","dashed"))


### Support Right & Logpopulation ####
zinb <- zeroinfl(executed_left ~ support_right + frontline + offset(logpopulation) 
                 + cntaffiliation + ugtaffiliation + border + sea + altitude  
                 + catholiccenter | support_right + frontline + population + ugtaffiliation + border + sea + altitude, 
                 data=nonunion, dist ="negbin"); summary(zinb)

data <- nonunion[,3:dim(nonunion)[2]] 
data <- apply(data, 2, mean,na.rm=T) 
median.values <- c("frontline","border","sea","catholiccenter") 
data[median.values] <- apply(nonunion[,median.values],2, median,na.rm=T)

count.vars <- c("support_right", "frontline", "cntaffiliation", "ugtaffiliation", "border", "sea", "altitude", "catholiccenter") 
zero.vars <-c("support_right", "frontline", "population", "ugtaffiliation", "border", "sea", "altitude")
vary <- "support_right"
logpop <- mean(nonunion$logpopulation, na.rm=T)[1]

setx <- XZout_Function(data=data, count.vars=count.vars, zero.vars=zero.vars, vary=vary)

sim.result <- ExpPred_Function.ZINB.log(zinb, x.out=setx$x.out, z.out=setx$z.out, logpop=logpop, par.simN=10^5, predict.simN=10^5)

plot(complevel, apply(sim.result$predict, 2, mean), type="b", ylim=c(0,15), xlab="Right Support", ylab="Number of Executions")
#title( main="Figure 1. Predicted Number of Executions by Competition", font.main = 1, cex.main=1)
lines(complevel, apply(sim.result$predict, 2, quantile, .025), type="l", lty="dotted")
lines(complevel, apply(sim.result$predict, 2, quantile, .975), type="l", lty="dotted")

lines(complevel, apply(sim.result$expect, 2, quantile, .025), lty="dashed")
lines(complevel, apply(sim.result$expect, 2, quantile, .975), lty="dashed")
legend(x="topleft", c("Predicted Value CI","Expected Value CI"), lty=c("dotted","dashed"))

################# FIRST DIFFERENCES FOR UNION INTERACTION MODEL ################
zinb <- zeroinfl(executed_left ~ support_right + union + support_right:union + frontline 
                 + offset(logpopulation) + cntaffiliation + ugtaffiliation + border + sea + altitude + catholiccenter | support_right 
                 + union + support_right:union + frontline + population + ugtaffiliation + border + sea + altitude, 
                 data=original, dist = "negbin") ; summary(zinb)

zinb$coeff$count
count.vars <- c("support_right","union", "frontline","cntaffiliation","ugtaffiliation",
                "border","sea","altitude","catholiccenter") # add in "support_right:union" later
zinb$coeff$zero
zero.vars <- c("support_right","union","frontline", "population","ugtaffiliation",
               "border","sea","altitude") # add in "support_right:union" later
data <- apply(na.omit(original[c("support_right","union", "frontline","cntaffiliation","ugtaffiliation",
                                 "border","sea","altitude","catholiccenter","population")]),2,mean,na.rm=T)
data[c("frontline","border","sea","catholiccenter")] <- 0 # median values
data
x.out1 <- XZout_Function(data=data, count.vars=count.vars, zero.vars=zero.vars,vary="support_right")$x.out[1,]
x.out1["support_right"] <- quantile(original$support_right,1,na.rm=T)
x.out0 <- x.out1
x.out0["support_right"] <- quantile(original$support_right,0,na.rm=T)
x.out1 <- c(x.out1,x.out1["support_right"]) # bind in interaction
x.out0 <- c(x.out0,x.out0["support_right"]) # bind in interaction
x.out0;x.out1

z.out1 <- XZout_Function(data=data, count.vars=count.vars,zero.vars=zero.vars,vary="support_right")$z.out[1,]
z.out1["support_right"] <- quantile(original$support_right,1,na.rm=T)
z.out0 <- z.out1
z.out0["support_right"] <- quantile(original$support_right,0,na.rm=T)
z.out1 <- c(z.out1,z.out1["support_right"]) # bind in interaction
z.out0 <- c(z.out0,z.out0["support_right"]) # bind in interaction
z.out0;z.out1


### first differences manual calculation
object <- zinb
par.simN <- 10^5
logpop <- mean(balcells$logpop,na.rm=T)
logpop

par <- object$optim$par
par.varcov <- -solve(object$optim$hess) # set variance covariance matrixs
# simulated parameters
sim.par <- rmvnorm(n = par.simN, par, par.varcov)

# get expected values distribution
mu.expect0 <- exp(sim.par[,1:length(x.out0)] %*% x.out0+logpop)
psi.expect0 <- (1/(1+exp(sim.par[,(length(x.out0)+1):(length(x.out0)+
  length(z.out0))]%*% -z.out0)))
mu.expect1 <- exp(sim.par[,1:length(x.out1)] %*% x.out1+logpop)
psi.expect1 <- (1/(1+exp(sim.par[,(length(x.out1)+1):(length(x.out1)+
  length(z.out1))]%*% -z.out1)))
sim.fd <- (1-psi.expect1)*mu.expect1 - (1-psi.expect0)*mu.expect0
mean(sim.fd)
quantile(sim.fd,0.975)
quantile(sim.fd,0.025)


####################################################################
############### FIRST DIFFERENCES FOR THE SEGREGATED DATA #######

nb <- glm.nb(executed_left ~ support_right + frontline + offset(logpopulation) + cntaffiliation 
             + ugtaffiliation + border + sea + altitude + catholiccenter,
             data=union) ; summary(nb)

x.out1 <- setx(nb,support_right=quantile(union$support_right,na.rm=T,1))
x.out1 
x.out1$frontline = x.out1$border = x.out1$sea = x.out1$catholiccenter <- 0 # set to medians (0)
x.out1

x.out0 <- setx(nb,support_right=quantile(union$support_right,na.rm=T,0))
x.out0 
x.out0$frontline = x.out1$border = x.out1$sea = x.out1$catholiccenter <- 0 # set to medians (0)
x.out0

fd.nb <- ExpPred_Function.NB.p2.log.fd(model=nb, x.out1, x.out0, 
                                       logpop=mean(union$logpopulation,na.rm=T), par.simN=10^5)

mean(fd.nb)
quantile(fd.nb,0.025)
quantile(fd.nb,0.975)




######################################################################################
##################### APPENDIX : Imputing Missing Data ###############################

impute <- balcells
impute.drop <- c("competition", "support_right", "compabs")
impute <- impute[,!(names(impute) %in% impute.drop)]
balcells.mi <- amelia(impute, m=5, idvars=c("id", "name"), noms=c("frontline", "border", "sea", "bombdum", "bombkilled10"), 
                      logs=c("population", "cntaffiliation", "ugtaffiliation"), sqrts=c("altitude"))

for(i in 1:5){
  balcells.mi$imputations[[i]]$support_right <- c(NA)
  balcells.mi$imputations[[i]]$support_right <- (100-balcells.mi$imputations[[i]]$support_left)
  balcells.mi$imputations[[i]]$competition <- c(NA)
  balcells.mi$imputations[[i]]$competition <- 1-((balcells.mi$imputations[[i]]$support_left-balcells.mi$imputations[[i]]$support_right)/100)^2
  balcells.mi$imputations[[i]]$compabs <- c(NA)
  balcells.mi$imputations[[i]]$compabs <- 1-abs((balcells.mi$imputations[[i]]$support_left-balcells.mi$imputations[[i]]$support_right)/100)
  balcells.mi$imputations[[i]]$union <- as.integer(0)
  balcells.mi$imputations[[i]]$union[(balcells.mi$imputations[[i]]$cntaffiliation>0|balcells.mi$imputations[[i]]$ugtaffiliation>0)] <- 1  
}

# Clean up #
rm(impute, i, impute.drop)
